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Introduction 


The  purpose  of  this  effort  was  to  investigate  the  sensitivity  of  various  shape 
recognition  algorithms  to  the  effects  of  shape  library  sampling  density  and 
imaging  noise  as  applied  to  the  task  of  helicopter  classification  and  orien¬ 
tation  determination. 

Section  1  introduces  the  first  of  two  shape  methods,  Fourier  descriptors, 
whose  performance  is  evaluated.  This  section  discusses  the  basic  method 
as  well  as  the  algorithms  for  the  efficient  computation  and  normalization 
of  the  shape  method  feature  vector. 

Section  2  introduces  the  second  shape  method,  Walsh  points.  The  def¬ 
inition  of  this  shape  method  is  stated.  Then  the  method  for  computation 
and  normalization  are  discussed. 

Section  3  discusses  the  procedure  used  to  generate  the  multiple  views 
of  the  reference  shapes  for  inclusion  into  the  reference  shape  feature  vector 
libraries. 

Section  4  describes  the  two  sets  of  experiments  to  determine  the  perfor¬ 
mance  of  the  two  algorithms  for  helicopter  fuselage  type  classification  and 
orientation  determination.  The  goal  of  the  experiments  is  to  determine  the 
sensitivity  to  library  sampling  density  and  imaging  noise. 

Section  5  interprets  the  results  of  the  experiments  described  in  Section 
4. 

Section  6  consists  of  a  proposed  scheme  for  the  library  procedures,  seg¬ 
mentation,  and  detection  of  the  helicopter  tip  path  for  rotor  blade  orienta¬ 
tion  determination. 

Section  7  provides  the  development  of  the  transformations  for  deter¬ 
mining  the  “unknown”  helicopter  orientation  given  the  orientation  of  the 
imaging  system  for  a  non-earth  centered  imaging  platform. 

The  results  of  the  entire  study  are  summarized  in  the  Conclusion. 

The  first  appendices  provides  the  raw  classification  and  angle  error  in¬ 
formation  including  the  classification  confusion  matrices  for  all  the  experi¬ 
mental  conditions  examined  in  the  study. 

The  last  appendix  is  a  listing  of  the  equations  required  to  compute  the 
elements  of  the  observed  object  orientation  relative  to  the  camera  reference 
frame. 
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1  Fourier  Descriptors 

1.1  Introduction 


The  first  method  used  in  recognizing  the  helicopter  fuseluge  silhouettes 
is  based  on  the  well  known  method  of  Fourier  series  analysis  [3],  [4],  [6], 
[lOj,  [15],  [16],  [17].  Having  been  provided  with  the  silhouette  of  an  object, 
the  contour  or  boundary  completely  specifies  the  tw'o-dimensional  shape. 
The  contour  can  be  parameterized  as  a  function  of  time  by  tracing  around 
the  boundary  in  counterclockwise  direction.  As  the  tracing  continues  the 
function  begins  to  repeat.  Since  this  function  is  periodic,  it  can  be  expanded 
into  a  Fourier  series.  Each  basis  function  of  the  complex  exponential  Fourier 
series  is  non-zero  almost  everywhere  over  each  period.  This  is  a  global 
method  of  shape  analysis.  Therefore,  when  segmentation  errors  become 
increasing  large,  the  recognition  accuracy  degrades. 

The  boundary  function,  q,  maps  the  real  number  line  into  the  complex 
plane.  The  projection  on  the  real  axis  is  the  x  component  and  the  orthog¬ 
onal  projection  on  the  imaginary  axis  is  the  y  component. 


MO  =  MO  +ty(t)  C  C. 

The  velocity  of  the  tracing  is  v(<)  =  So,  the  speed  of  tracing  is 

equal  to  the  magnitude  of  the  velocity,  i.e. 

■  ,  ,(f)  12  dK(t)dy(t) 

— IT- 


If  we  assume  that  the  tracing  takes  place  at  a  constant  speed  ( uniform 
tracing),  then 


v 


dl(t) 

dt 


1/2 

=  constant. 


Since  the  speed  is  constant  and  the  period  of  the  trace  is  T,  then  T  = 
vL,  where  L  is  the  total  arc  length  once  around  the  contour.  Assuming 
q(/)  is  a  continuous,  bounded,  periodic  function  with  finite  arc  length  (i.e. 
rectifiable),  -7(f)  can  be  expanded  in  the  Fourier  series 


MO 


n  -  —  oc 


2 


where 


1.2  Calculation 

Since  the  images  are  discrete,  connecting  the  contours  of  the  pixels  around 
the  boundary  produces  a  polygon.  The  Fourier  coefficients  can  be  com¬ 
puted  directly  from  the  increments  in  time  (arc  length)  and  position  as  the 
polygon  is  traced.  This  is  called  the  direct  Fourier  transform  (DFT).  The 
following  derivation  follows  along  the  lines  of  [6],  [lOj. 

Let  tiie  increment  in  position  be  the  complex  number  A 7,,  and  the 
increment  in  time  will  be  A  tt  =  |  A  7,  |  =  yArr,2  +  Ay,2.  Since  7  is  piecewise 
linear,  its  derivative,  7,  is  piecewise  constant.  So,  let 

i{t)  =  £ 

n=  -  00 


be  the  Fourier  series  for  7.  But  the  series  for  7  can  be  differentiated  term 
by  term  so  that 

+oc  .2  itn 


So, 

Now, 


7 


„  .2tt 

Pn  —  1  rp  Cn 


p  ly  _  e-*¥‘r-)  (t  -1_) 

T  &tpK  1  1  27m' 


7=  £ 


■  "7  ( 


—  ft.,  n  #  0. 

27rn 


.  T 


We  can  then  write 


w 


here 


47r2n2  ~  Afp 

p 

tv  -  £  AL,  -  0, 


e  ,2-r ,  n  ^  0, 


(1) 


1 


3 


a,  : 

0 

1  2 

3 

4 

5  6 

7 

A  ^  : 

1 

1  +  i  i 

-1  +  i 

-1 

-1-i  -i 

1-i 

At  : 

1 

y/2  1 

V~2 

1 

y/2  1 

V2 

Table  1:  Contour  line  increments  for  the  8-direction  chain  code. 

K  is  the  number  of  sides  on  the  polygon,  and  T  =  tf (  —  period.  For  the  c0 
term, 

1  * 

co  —  T  /7p-i)A£p.  (2) 

1  p= i 

For  the  discrete  images,  the  contours  are  represented  by  the  eight- 
direction  chain  code.  So,  the  A ^  and  At,  can  be  determined  using  a 
lookup  table  indexed  by  the  chain  code. 

The  computational  complexity  is  of  order  NK,  where  N  is  the  number 
of  coefficients  calculated  and  K  is  the  number  of  links  in  the  chain  coded 
contour. 

1.3  Normalization 

It  is  usually  desirable  to  compare  shapes  independent  of  size,  orientation, 
and  starting  point.  To  compare  two  shape  boundaries  using  their  Fourier 
descriptors,  it  is  necessary  to  scale,  rotate,  and  shift  their  shapes  in  order 
to  allow  the  “best”  fit  possible.  This  operation  normalizes  the  Fourier 
descriptors  for  the  unknown  shape  to  an  optimum  orientation.  It  has  been 
shown  [15]  that  the  optimum  scale,  rotation  angle,  and  relative  starting 
point  shift  can  be  obtained  to  minimize  the  mean-squared  error  as  the 
criterion. 

To  reduce  the  computational  burden  to  normalize  the  coefficients,  many 
others  have  suggested  that  one  of  several  suboptimum  methods  be  used. 
Instead  of  normalizing  the  coefficients  of  the  unknown  differently  for  each 
template,  the  coefficients  are  normalized  to  a  “standard”  orientation  inde¬ 
pendent  of  the  template. 

First  the  object  is  translated  to  the  origin  by  setting  Co  =  0.  Because  the 
fundamental  or  C\  has  always  been  observed  to  be  the  largest  (n  yt  o),  all 
the  coefficients  are  scaled  by  dividing  by  |  C\  |.  Next,  the  shape  is  rotated 
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and  the  starting  point  shifted.  Most  of  the  variations  among  the  suboptimal 
methods  described  in  the  literature  occur  in  how  this  rotation  and  starting 
point  shift  is  obtained  [12],  [  1 3] ,  [I6i,  [  17]. 

The  following  algorithm  is  used  in  this  study  to  normalize  the  shape 
coefficients: 

1.  Set  cq  =  0. 

2.  Scale  so  that  I  c\  |=  1. 

3.  Find  the  coefficient  ck  that  is  next  largest  (j  k  —  1  j<  5,  k  ^  0,  1.)  If 
j  k  —  1  |>  5,  then  let  k  —  2. 

4.  Rotate  and  shift  the  starting  point  so  that  lc\  =  0  and  Zc*  =  0,  i.e. 


Cn 


Cn  tfnCi  4-a) 

I  c  ’ 

C1  I 


where 


/Ci  -  /ck 


{k/ci  -  Z ck ) 
k  -  1 


The  object  -7 (t)  =  Cie'^‘  f  cke ,l^t  has  ]  k  —  1  [-fold  symmetry.  So,  there 
are  |  k  —  1  ]  rotations  and  relative  starting  point  shifts,  multiples  of 
that  will  satisfy  the  zero  phase  condition.  So,  rotate  and  shift  the  starting 
point  so  that 

c„(m)  -  c'„em‘ri  _‘rrr,  m  —  0,  A:  —  1  |  — 2). 


to  maximize  the  following  function: 

|  5?{cn(m)}  |  . 

n 

This  criterion  effectively  chooses  the  normalization  that  orients  the  contour 
so  that  the  axis  of  one  of  the  main  lobes  of  the  |  k  —  1  j-fold  fundamental 
shape  c  1  e,£  +  cke'kt  is  along  the  positive  x-axis  and  the  starting  point  on  the 
contour  corresponding  to  that  lobe  where  the  contour  is  farthest  from  the 
origin.  In  order  to  reduce  the  number  of  rotations  to  perforn  and  calculate 
the  criterion  only,  j  k  —  I  j<  5  is  allowed.  If  |  k  —  1  |>  5,  then  k  —  —  1. 

The  correct  recognition  of  a  shape  is  very  sensitive  to  this  rotation  and 
starting  point  normalization.  In  order  to  improve  the  classification  ac¬ 
curacy,  multiple  sets  of  Fourier  coefficients  can  be  used  in  classifying  the 
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shape.  The  descriptor  for  me  best  normalization  above  is  used.  In  addition, 
if  one  of  the  multiple  rotations  has  a  value  of  its  criterion  that  is  within 
95%  of  the  best,  this  normalization  is  also  used.  To  reduce  the  possibilities 
that  the  wrong  coefficient  is  used  in  the  normalization,  the  third  largest 
coefficient  is  also  used  if  it’s  magnitude  is  within  95%  of  the  second  largest 
coefficient.  Again,  the  next  best  normalization  based  on  this  third  largest 
coefficient  is  also  used  if  the  criterion  is  95%  of  the  best  normalization 
based  on  this  new  coefficient.  So  in  all  there  can  be  as  many  as  four  sets  of 
coefficients  for  each  unknown  shape.  For  the  library  features  only  the  best 
normalization  is  used.  One  of  the  four  possible  normalizations  is  likely  t~ 
match  the  proper  template  even  when  some  noise  is  present. 

1.4  Fourier  Descriptor  Feature  Vector 

The  feature  vector  was  formed  by  calculating  and  normalizing  the  coeffi¬ 
cients.  After  normalization,  the  Cq  and  components  did  not  carry  any 
shape  information  and  were  dropped  from  the  feature  vector.  The  fea¬ 
ture  vector  was  then  formed  by  listing  the  real  and  imaginary  parts  of  the 
remaining  components  as 

/  =  («{c- 1},  5{c_ 

3?{c-2},5{c_2},. .  .,9?{cN/2},Q,{ciV/2})T. 

This  is  the  full  feature  vector.  The  32  coefficients,  c_i6  to  c  15,  are  used  for 
the  experiments. 

2  Walsh  Points 

2.1  Introduction 

The  Walsh  functions  have  often  been  used  as  the  basis  for  functional  approx¬ 
imation  because  the  piecewLe  constant  form  leads  to  efficient  computation 

3  .  Also,  since  the  basis  functions  are  non-zero  over  the  entire  interval 
of  defmition,  this  is  a  complete  global  shape  analysis  method.  Instead  of 
the  transform  coefficients,  this  method  uses  the  basis  functions  to  formu¬ 
late  the  computation  of  Walsh  points  for  the  shape  feature  set.  1  he  Wai.  n 
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points  are  a  collection  of  points,  xk  +  iyk,  k  =  0, 1, 2, . . . ,  2m_1 ,  that  have 
a  property  derived  from  the  fact  that  the  Walsh  functions  were  used  in 
formulating  their  calculation.  The  projections  x(t)  and  y(t)  are  approxi¬ 
mated  to  within  a  prespecified  degree  of  accuracy  by  a  piecewise  constant 
function,  a  truncated  Walsh  series. 

The  Walsh  functions,  H ',(<),  are  products  of  the  Rademacher  functions, 
r,(t).  The  Walsh  functions  for  a  complete  orthogonal  basis  set  complete 
among  square  integrable  functions  on  the  interval  0,  1 1 .  The  Walsh  func¬ 
tions  are  defined  as 

H'o(0  -  1. 

IMO  rnl  +  l{t)rri2-i  l  '  ■  •  •  '  +  1  (0  ' 

where  n  >  1  is  expressed  in  binary  as 

n  =  2n‘  +  2n-  -*•  . . .  +  2np 
and  the  integers  n,  are  ordered  such  that 


rii  <  n2  <  .  . .  <  T2p. 


The  Rademacher  functions  are 


ro(0 

ri(0 

ri{t  +  l) 

rki  1 (0 


_  |  +M  t  ;0,  D 

=  n(0 

=  ri(2kt),k  —  0, 1, 2, . . . 


In  the  following  only  the  x  projection  of  the  boundary  calculation  will  be 
discussed.  The  same  results  follow  for  the  y  projection. 

The  Walsh  series  for  x(t)  is 

OO 

x(<)  = 

n~  0 
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where 

a,  =  [  x(t)WAt/T)dt  =  T  [  x(Tt)WAt)dt 

Jo  Jo 

T  =  {total  arc  length  along  7(f)  =  x{t)  4-  iy(<)}. 

Each  a,  is  made  up  of  the  area  under  the  portions  of  x(t)  added  and  sub¬ 
tracted  together.  In  approximating  x(t),  only  a  finite  number  of  terms  in 
the  Walsh  series  are  retained  -  say  the  first  2m.  If  the  series  is  truncated 
to  a  finite  number  of  terms,  there  is  an  approximation  error.  It  would  be 
expected  that  as  N  =  2m  increases,  the  closer  x2m(£)  approximates  x(t). 
Let 

2m  —  1 

x2 ™(0  -  £  akWk{t/T). 

k= o 

Since  the  function  z2m  ( t )  (or  y2m  (f ))  describes  a  piecewise  constant  function, 
there  are  2m  discrete  points,  xk  +  iyk,  that  are  the  sum  of  the  heights  of 
the  Walsh  functions  over  each  interval,  T ,  ^T),  k  =  0,1,...,  2m  —  1, 
which  approximate  the  x  (or  y)  projection  of  7(f)  =  x(f)  +  iy[t). 

It  was  mentioned  earlier  that  the  Walsh  functions  allow  for  an  efficient 
computation  of  the  function  x2m(t)  (or  y2m(().)  The  ith  Walsh  function  can 
be  written  as 

*=0 

where 

-»*  =  »’■(<)  =  ±i. 

The  matrix  of  7,*,  H  =  ( //,>]  =  [7^],  is  the  Hadamard  matrix  with  2m  rows 
and  columns. 

The  ith  coefficient  of  the  Walsh  series  is 

=  T  fQ  x(Tt)\Vi[t)dt  =  T  YL  Aik  Jq  ziT^x^^dt 

2m  —  1  II 

=  X*  Aik  I kr  x(t)dt 

k=o  J  r" 

2m  —  1 

o.  —  'y  ^  7i*^it)  (3) 

k  =  0 
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where 


,  IHiit 

T  Jl  X(Tt)dt=Jkr  xit)dL 


*T 
2  m 


Ax  — 


ylxjt  is  the  area  under  x[t)  over  the  basic  interval  ^k^JT).  By  letting 


/  AQz  ^ 

-^lx 

V  Av" -ix  y 

f  ao  \ 

«1 


be  the  Zrea  vector  and 


&r  — 


V  &2"*-l  / 

be  the  sequericy  vector,  we  can  write 


d*,  —  H  At. 


The  Truncated  Walsh  series  is 


*i~(0  = 

=  I>.(l :‘-r..xl*.^l«/r) 


*=o  V  .=o  ' 

Yl  c*^!-24t.-2w  )(t/T)' 


k  *4J 
2 m  ■  2,ri 


M/T) 


where 


2m  -  1 

i  — 0 


(4) 
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The  constants  c*  can  be  found  by  first  forming  the  vector  xfm  as 


£2" 


Co 
Cl 

V  C2">-1  J 


and  noticing  that  as  in  equation  (1)  that  equation  (2)  can  be  written  as 

X  2  m  =  II  ^  Cl  z 

where  HT  is  the  transpose  of  the  Hadamard  matrix,  H . 

Since  II  is  symmetric  II  =  HT  and 

X  2 m  -  H2Ax. 


But  because  the  Walsh  functions  are  orthogonal,  it  follows  that  H 2  =  2m. 
Hence, 

X2m  ~  2m  Az. 

So,  the  truncated  Walsh  series  gives  a  piecewise  constant  function  whose 
heights  are  proportional  to  the  area  under  the  projection  on  each  basic 
interval,  i.e. 

2m  - 1 

x2™(0  =  2m  y:  AkxXikr  jk+pT* (Q. 

I  2m  I  2m  1 


2.2  Normalization 

In  order  to  compare  the  Walsh  points  feature  vector  of  an  unknown  to  a 
library  prototype,  normalization  with  respect  to  translation,  scale,  rotation, 
and  shift  in  starting  point  must  be  accomplished. 

The  traced  contour  from  the  image  grid  is  stored  using  Freeman  chain 
code  links.  When  the  shape  features  are  to  be  extracted,  the  chain  code  is 
converted  to  a  complex  vector  7  =  £+  iy.  The  increments  in  the  arc  length 
for  the  piecewise  linear  segments  connecting  the  grid  points  are  calculated. 
The  c0,  Ci,  and  c_i  Fourier  coefficients  of  the  boundary  are  calculated  to 
provide  the  global  information  necessary  for  translation,  scale,  and  rotation. 
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The  translation  is  —  c0.  The  scale  factor  5  =|  c\  |.  The  rotation  angle  is 

(  7o  ^ 

,  .  The  complex  vector  -7  = 


the  ith  component  is 


is  normalized  to  7,  where 


V  7k  7 


7,  =  —  (nr,-  —  Co)^,  and  T  =—. 


This  normalization  moves  the  center  of  the  shape  to  the  origin.  It  scales  and 
rotates  the  shape  so  that  the  major  axis  corresponding  to  the  fundamental 
ellipse  is  along  the  x-axis  and  has  a  length  of  2.  The  starting  point  is  moved 
by  finding  the  kth  point,  -7*,  which  is  closest  to  one  of  the  two  points  where 
the  fundamental  ellipse  crosses  the  x-axis.  An  additional  180c  rotation  is 
applied  if  necessary  to  place  -7*  in  the  right  half  plane.  So,  the  normalized 
Walsh  points  form  the  vector 


7  = 


Mk  \ 

7  k -1 
7o 

V  7k- 1  J 


shift  7 
7- 


or  5  =  7  •  e* 2  ■  Then  the  Walsh  points  are  calculated  by  computing 

K 

the  Area  vector  and  multiplying  by  2m. 


2.3  Walsh  Points  Feature  Vector 

The  feature  vector  used  in  the  experiments  was  formed  by  simply  listing 
the  Fourier  points  real  and  imaginary  parts  in  order,  i.e. 

/  =  (^oi  yiu  xi,  yi  i  5  xk-  1 1  vk-  1  )r 

For  the  experiments  the  full  feature  vector  consisted  of  32  Walsh  points. 
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3  Shape  Library  Generation 

The  focus  of  this  study  is  to  determine  the  sensitivity  of  the  shape  recog¬ 
nition  algorithms  to  the  density  of  library  entries  and  to  imaging  noise. 
The  algorithms  are  examined  for  their  performance  for  correct  helicopter 
fuselage  type  recognition  and  orientation  determination. 

A  list  of  shape  features  are  extracted  from  the  outer  contour  of  the  sil¬ 
houette  of  a  helicopter  fuselage.  These  features  are  derived  from  computer 
models  of  the  helicopter  fuselage.  A  computer  program  is  used  to  obtain  a 
three-dimensional  model  of  the  solid  object.  Additional  computer  graphics 
programs  are  used  to  rotate  the  internal  model  and  produce  a  computer 
generated  image  of  the  solid  object.  Therefore,  for  each  shape  that  the 
algorithms  are  to  recognize,  it  is  necessary  to  enter  into  the  computer  its 
three-dimensional  description.  A  three-dimension  solid  modelling  program, 
Mac3D  [11],  generates  the  solid  models.  This  program  allows  the  user  to 
create  a  three-dimensional  model  by  specifying  the  size  and  orientation 
of  simple  three-dimensional  objects.  These  objects  include  spheres,  cones, 
parallel-pipeds,  tori,  and  surfaces  of  revolution.  Very  complex  objects  can 
be  modelled  by  combining  these  simple  shapes.  A  good  approximation  to 
the  real  world  object  shape  can  be  obtained  in  this  manner. 

The  helicopter  rotors  blades  are  not  included  in  the  modelling  of  the  he¬ 
licopter  fuselage.  This  is  because  experiments  in  segmenting  the  helicopter 
fuselage  have  indicated  that  the  rotor  blades  are  in  fact  difficult  to  extract 
from  the  image.  This  topic  is  discussed  in  detail  in  Section  7. 

The  three-dimensional  specifications  for  three  military  helicopters  are 
entered  into  the  solid  modelling  program  for  this  study.  The  descriptions 
are  taken  from  Jane’s  All  the  World’s  Aircraft  [9].  The  helicopters  choosen 
are:  the  Bell  Model  205  UH-1  Iroquis,  the  McDonnell  Douglas  500MD 
Defender,  and  the  Sikorsky  UH-60A  Black  Hawk. 

Once  the  three-dimensional  object  description  is  available,  the  image  of 
the  helicopter  in  any  arbitrary  orientation  can  be  obtained.  The  object  ro¬ 
tation  transformations  are  discussed  in  Section  X.  Before  the  classification 
and  orientation  experiments  are  performed,  a  library  of  views  is  generated 
for  each  helicopter  type.  For  each  library  for  the  three  types  of  helicopters, 
the  following  procedure  is  carried  out.  For  each  view  to  be  stored  in  the 
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library,  the  object  is  rotated  by  the  specified  rotation  about  the  x  and  y 
axis,#*  and  0V,  respectively.  The  object  coordinate  points  are  then  projected 
onto  the  image  plane.  The  image  formed  is  the  silhouette  of  the  helicopter 
fuselage  viewed  in  the  specified  orientation.  The  the  contour  of  this  sil¬ 
houette  is  then  traced  using  a  8-direction  chain-code.  The  chain-codes  for 
each  stored  view  are  collected  into  the  library  for  that  helicopter  type.  The 
angle  orientation  information  is  retained  in  a  special  data  record  with  each 
view  for  latter  use.  The  chain-code  libraries  are  then  processed  using  one 
of  the  two  shape  recognition  algorithms  to  extract  the  shape  features.  The 
normalized  shape  features  for  each  helicopter  type  are  stored  together  in  a 
feature  vector  library. 

Once  the  libraries  containing  the  shape  features  for  each  helicopter  have 
been  generated,  experiments  can  be  performed.  The  performance  of  each 
helicopter  shape  is  determined  by  classifying  labelled  ”unknown”s  from 
each  helicopter  type  viewed  at  random  orientations.  Shape  feature  vectors 
are  classified  using  a  nearest  neighbor  template  matching.  The  normal¬ 
ized  shape  feature  vector  is  compared  to  every  feature  vector  in  each  of 
the  three  libraries.  The  unknown  is  assigned  to  the  helicopter  type  and 
orientation  corresponding  to  the  closest  match.  The  distance  measure  used 
for  the  template  matching  is  the  sum  of  the  squared  differences  between 
corresponding  feature  vector  elements,  i.e. 


d 


h 


i=i 


where  /,  1  is  the  ith  element  in  the  jth  feature  vector  of  the  Lth  helicopter 
normalized  shape  feature  vector  and  /,u  is  the  ith  element  of  the  unknown 
normalized  shape  feature  vector. 

Libraries  for  each  helicopter  type  are  generated  at  various  angle  sampling 
densities.  The  views  that  make  up  each  library  are  taken  from  equally 
spaced  values  in  the  x  and  y  orientation  angles.  The  x  and  y  angle  rotations 
range  from  0  to  -n  and  —  |  to  |,  respectively.  For  each  helicopter  type, 
four  different  libraries  are  generated  for  the  purpose  of  determining  the 
sensitivity  of  classification  and  orientation  accuracy  to  the  density  of  the 
library  .  The  four  libraries  for  each  helicopter  type  contain  either  5x3,  7x5. 
9x7,  or  11x9  views. 
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4  Experimental  Procedure 

Two  groups  of  experiments  are  performed  to  determine  the  performance 
of  the  two  shape  recognition  algorithms  as  measured  by  helicopter  fuselage 
type  classification  and  orientation  angle  accuracy.  The  purpose  of  the  first 
group  experiment  is  to  determine  the  sensitivity  of  the  algorithm  perfor¬ 
mance  as  a  function  of  the  density  of  reference  library  entries.  The  second 
group  of  experiments  is  to  determine  the  sensitivity  to  imaging  noise. 

For  the  first  group  of  experiments,  libraries  for  each  helicopter  type  are 
generated  at  various  angle  sampling  densities.  The  views  that  make  up  each 
library  are  taken  from  equally  spaced  values  in  the  x  and  y  orientation 
angles.  The  x  and  y  angle  rotations  range  from  0  to  — 7r  and  —  *  to  * , 
respectively.  For  each  helicopter  type,  four  different  libraries  are  generated. 
The  four  libraries  for  each  helicopter  type  contain  either  5x3,  7x5,  9x7,  or 
11x9  views  corresponding  to  orientation  angle  sampling  intervals  of  30°  v 
60%  25.7°  x  36°,  20°  x  25.7°,  and  16.4°  x  20°.  The  contours  for  each  of 
the  three  helicopter  types  at  each  of  the  four  library  sampling  densities  are 
shown  in  Figures  1-12.  In  the  libraries  shown  in  the  figures,  the  variation 
in  0V  are  horizontal  and  the  variations  in  dz  are  vertical.  The  (0r,0J  -  (0,0) 
is  at  the  top  center. 

For  each  classification  experiment,  the  unknowns  are  compared  to  a  li¬ 
brary  for  each  helicopter  type,  all  at  the  same  library  density.  In  each 
experiment,  03  random  views  of  each  of  the  three  helicopter  types  are  clas¬ 
sified.  The  classification  accuracy  and  statistics  on  the  angle  orientation 
errors  are  computed  and  tabulated.  This  procedure  is  carried  out  using  for 
both  the  Fourier  descriptors  and  the  Walsh  points.  A  small  sample  of  the 
“unknowns"  arc  shown  in  Figures  14  -  16. 

The  purpose  of  the  second  set  of  experiments  is  to  examine  the  efTect 
of  imaging  noise  on  classification  accuracy  and  orientation  determination 
of  a  helicopter  fuselage.  In  addition,  the  experiments  are  performed  with 
the  feature  vector  shape  libraries  at  each  of  the  four  sampling  densities.  As 
before,  the  outer  contours  for  63  random  views  of  each  helicopter  type  are 
generated  using  the  solid  modelling  program  followed  by  object  coordinate 
rotation  and  projection.  The  outer  contour  of  the  rotated  and  projected 
.vurfac edges  of  the  solid  object  are  traced.  Next,  the  interior  of  the  contour 
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Figure  2:  McDonnell  Douglas  500MD  Defender  library  of  5x3  (15)  views. 
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Figure  3:  Sikorsky  UIIGOA  Black  Hawk  library  of  5x3  (15)  views. 
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Figure  5:  McDonnell  Douglas  500MD  Defender  library  of  7 \5  (35)  views. 
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Figure  6:  Sikorsky  UH60A  Black  Ilawk  library  of  7x5  (35)  views. 
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Figure  7:  Bell  Model  205  I  II- 1  Iroquois  library  of  9x7  (63)  views. 
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Figure  8:  McDonnell  Douglas  500MD  Defender  library  of  9x7  (63)  views. 
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Figure'  9:  Sikorsky  UIIGOA  Black  Hawk  library  ot  9x7  (03)  views. 
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Figure  10:  Bell  Model  205  UH-1  Iroquois  library  of  11x9  (99)  views. 
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Figure  11:  McDonnell  Douglas  500MI)  Defender  library  of  11x9  (99)  views 
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Figure  12:  Sikorsky  UH60A  Black  Hawk  library  of  11x9  (99)  views. 
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Figure  13:  Example  contours  for  1  o  “unknowns  -  Hell  Model  20 


Figure  14:  Example  contours  for  15  “unknows”  -  McDonnell  Douglas 
500MD. 
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Figure  15:  Example  contours  for  15  “unknowns"  -  Sikorsky  UHGOA. 
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is  filled  with  a  level  of  191  and  the  outside  area  (background)  the  value 
61.  White  Gaussian  noise  is  then  added  to  simulate  image  noise.  The 
image  is  then  thresholded  at  a  level  of  110.  The  outer  contours  in  the 
image  are  then  traced  and  the  longest  contour  retained  to  describe  the 
noisy,  helicopter  fuselage  silhouette,  outer  contour.  Two  additional  sets 
of  labeled  “unknown”  helicopter  contours  are  generated  in  this  manner. 
Each  set  corresponds  to  a  object  intensity  level  above  background  to  noise 
variance  signal-to-noise  ratio,  i.e. 

SNR  —  20log(  — ), 
a 

where  A  =  (object  gray  level  intensity  -  background  grey-level  intensity) 
and  a  —  the  standard  deviation  of  the  noise  process.  The  classification 
and  orientation  accuracy  experiments  are  then  performed  using  unknowns 
at  the  three  signal-to-noise  ratios,  infinity,  lOdB,  and  3  dB  corresponding 
to  noise  a  of  0,  20.238,  and  45.310,  respectively.  Examples  of  the  images 
and  contours  at  each  signal-to-noise  ratio  for  the  Bell  N1D500  helicopter  are 
shown  in  Figures  16  -  21. 

5  Recognition  Accuracy  and  Orientation  De¬ 
termination  Results 

The  ability  of  the  shape  recognition  algorithms  to  assign  the  “unknown” 
shape  contour  to  the  correct  helicopter  fuselage  type  are  tabulated  in  Ta¬ 
ble  2  and  shown  in  Figure  22.  The  plot  and  table  show  the  classification 
accuracy  under  all  the  tested  experimental  conditions  including  variation 
in  library  density  and  image  signal-to-noise  ratio.  Under  all  conditions 
the  Fourier  descriptor  shape  recognition  algorithm  out  performs  the  Walsh 
points  recognition  algorithm.  The  margin  is  always  at  least  and  additional 
12%.  The  performance  of  both  algorithms  degrade  significantly  when  the 
number  of  library  entries  falls  below  35  (7x5).  With  a  moderate  decrease 
in  the  signal-to-noise  ratio,  little  degradation  in  classification  accuracy  oc¬ 
curs.  However,  when  the  signal-to-noise  ratio  is  decreased  to  only  3dB,  the 
Fourier  descriptor  accuracy  drops  some  8%.  The  Walsh  points  algorithm 
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Figure  18:  He!!  M  1)500  silhouette  a:u 


Figure  19:  Bell  MD500  contour  with  no  imaging  noise. 
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Figure  20:  Bell  MD500  noisy  contour  at  SNR  =  lOdB. 


is  more  sensitive  and  drops  by  at  least  13%  when  the  signal-to-noise  ratio 
decreases  from  10  dB  to  3  dB.  The  sensitivity  curves  indicate  a  non-linear 
effect.  The  orientation  normalization  used  in  both  algorithms  is  a  non¬ 
linear  decision  process.  This  would  explain  this  non-linear  phenomenon. 
At  a  signal-to-noise  ratio  of  3  dB,  the  effects  of  partial  contours  due  to 
segmentation  errors  also  just  begins  to  occur.  Under  such  conditions,  both 
the  performance  of  Fourier  and  Walsh  methods  deteriorates  rapidly  }6.. 

Figures  23  thru  2G  and  Tables  3  and  4  show  the  sensitivity  of  the  al¬ 
gorithms  to  the  experimental  conditions  as  measured  by  orientation  angle 
error.  Plots  for  the  mean  and  median  x  and  y  angle  errors  are  shown. 
Again,  under  almost  all  conditions,  the  Fourier  descriptor  method  out  per¬ 
forms  the  Walsh  points  method.  The  angle  errors  for  the  Walsh  method 
being  approximately  twice  as  large  as  those  for  the  Fourier  method.  For  the 
Fourier  descriptors,  the  y  angle  error  is  approximately  one-half  the  library 
y  angle  sampling  interval.  This  is  the  best  that  is  likely  to  be  obtained  for 
any  shape  method  without  employing  some  interpolation  between  library 
views.  For  the  Walsh  points  method,  the  y  angle  errors  are  twice  that  of 
the  Fourier  method. 

For  the  Fourier  descriptor  method,  the  mean  errors  for  the  x  angle  ori¬ 
entation  determination  is  approximately  two  and  a  half  times  as  large  as 
the  library  x  angle  sampling  interval.  The  mean  x  angle  error  for  the  Walsh 
points  method  is  as  much  as  four  times  the  library  sampling  interval.  The 
cause  of  the  large  mean  x  angle  errors  is  first  indicated  by  examining  the 
median  x  angle  errors,  f'or  example,  the  mean  x  angle  error  for  the  Fourier 
descriptor  method  at  the  11x9  library  density  is  25c,  but  has  a  median  error 
of  only  8°.  Inspecting  the  few  individual  cases  where  there  is  a  large  x  angle 
error,  shows  that  the  same  fuselage  at  a  displacement  angle  of  90”  is  being 
incorrectly  assigned.  The  number  of  such  instances  is  few  but  the  90r  dis¬ 
crepancy  is  large  enough  to  greatly  bias  the  sample  mean  error.  The  effect 
of  such  errors  could  be  mitigated  in  an  operating  scenario,  by  smoothing 
the  angle  orientations  after  locking  on  to  the  target  trajectory. 

The  median  angles  for  the  Walsh  points  method  are  also  small  indicating 
the  same  problem.  However,  the  performance  is  already  inferior  to  that  of 
the  Fourier  method. 

In  summary,  the  overall  performance  of  the  Fourier  descriptor  method 
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is  superior  to  that  of  the  Walsh  points  method.  Both  methods,  however, 
exhibit  large  mean  x  angle  errors.  However,  most  of  the  error  is  due  to  a  few 
large  miss  assignments  in  orientation  that  can  be  post-filtered  by  a  tracking 
algorithm  which  incorporating  the  limits  of  the  helicopter  dynamics. 


6  Rotor  Orientation  Determination 

In  this  section,  a  scheme  for  segmenting  and  subsequently  determining  the 
helicopter  rotor  orientation  relative  to  the  helicopter  fuselage  will  be  dis¬ 
cussed.  If  rotor  orientation  can  be  determined  from  the  imagery  in  a  timely 
fashion,  it  should  be  possible  to  provide  a  more  accurate  prediction  of  the 
future  trajectory  of  the  helicopter. 

In  pursuit  of  the  above  goal,  video  imagery  of  a  maneuvering  helicopter 
was  obtained  and  digitized.  Several  algorithms  wrere  applied  in  an  attempt 
to  segment  the  rotor  blade  “disk”  from  the  image.  All  the  methods  failed  to 
produce  a  useable  segmentation.  It  is  this  same  fact,  however,  that  makes 
it  possible  to  apply  the  global  shape  methods  in  Sections  1  and  2  of  this 
report.  The  difficulty  in  segmenting  the  rotor  blades  reduces  the  confusion 
in  providing  a  complete  shape  of  the  helicopter  silhouette.  However,  further 
investigation  in  to  the  acquisition  of  the  data  on  the  video  tape  indicates 
that  an  alternative  should  exist. 

The  data  recorded  on  the  video  tape  was  acquired  using  a  vidicon  cam¬ 
era.  This  sensor  type  experiences  image  lag.  Therefore,  the  image  of  the 
individual  rotor  blades  were  integrated  temporally  resulting  in  a  single  very 
low  contrast  elliptical  disk  of  the  rotor  blade  tip  path.  (For  most  orienta¬ 
tions  of  the  helicopter  fuselage.)  By  using  a  different  sensor  type,  e.g.  a 
CCD  camera,  the  effects  of  temporal  averaging  can  be  significantly  reduced. 
Such  a  sensor  operating  at  standard  rate  of  60  fields  per  second  would  image 
the  rotor  blade  through  only  0.1  of  a  revolution.  (Assuming  rotor  radius  of 
27  feet  and  a  tip  speed  of  Mach  1.)  Having  a  short  integration  time  requires 
that  the  camera  iris  be  open  as  wide  as  possible,  thereby  reducing  the  depth 
of  focus.  This  forces  the  constraint  that  the  camera  already  be  locked  on 
the  target  in  a  tracking  mode.  Subsequently,  the  contrast  between  nominal 
background  and  the  rotor  blades  should  be  significantly  enhanced.  Making 
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Figure  22:  Performance:  Classification  Accuracy 
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Figure  23:  Performance:  Mean  X  Angle  Error 
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Figure  24:  Performance:  Mean  Y  Angle  Error 


Table  4:  Performance:  Median  Angle  Error 
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Figure  25:  Performance:  Median  X  Angle  Error 


these  assumptions,  the  following  procedure  is  proposed  to  determine  rotor 
orientation. 

First,  the  shape  libraries  for  each  hehcopter  type  is  enhanced  to  contain 
with  each  view  the  x  and  y  displacement  to  the  rotor  center  of  rotation.  In 
addition,  with  each  library  is  included  the  eight  vertices  of  a  rectangular 
prism  which  is  a  bounding  surface  for  the  envelope  of  the  rotor  blades. 
These  coordinates  and  displacements  are  scaled  to  the  reference  image  size 
used  in  the  generation  of  the  shape  libraries.  With  each  library  entry  is 
the  scale  required  for  the  shape  vector  normalization.  The  ratio  between 
the  tracked  “unknown”  contour  normalization  scale  and  the  stored  scale  for 
the  assigned  library  entry  orientation  provides  the  scale  factor  required  to 
scale  the  displacement  vector  and  bounding  surface  vertices  for  the  actual 
scale  of  the  imaged  target.  The  rotation  orientation  angles  6Z,  6 fc,  and  9Z 
are  then  used  to  rotate  the  bounding  surface  coordinates  and  then  project 
them  onto  the  image  plane.  The  rotor  center  of  rotation  coordinates  are 
rotated  in  the  image  plane  by  6Z.  The  outer  contour  of  the  scaled  rotated, 
displaced,  and  projected  bounding  surface  will  provide  a  limited  region  of 
interest  on  which  to  perform  the  rest  of  the  operations.  Limiting  the  region 
of  interest  will  reduce  both  the  total  processing  time  required  as  well  as 
limit  the  interfere  of  clutter  in  the  image. 

Once  the  limits  of  the  imaged  bounding  surface  is  obtained,  the  magni¬ 
tude  and  angle  of  the  gradient  is  computed  for  all  points  inside  the  limits, 
but  not  on  the  silhouette  of  the  helicopter  fuselage.  A  Hough  transform  [8’ 
is  then  computed  to  detect  the  lines  delineating  the  rotor  blades.  In  com¬ 
puting  the  Hough  transform  only  pixels  inside  the  bounding  region  and  not 
inside  the  contour  of  the  helicopter  fuselage  (which  is  readily  available)  are 
plotted  in  Hough  space.  This  is  to  reduce  the  inference  of  linear  features  in¬ 
terior  to  and  delineating  the  fuselage  contour.  The  lines  detected  in  Hough 
space  are  then  followed  inside  the  limits  of  the  bounding  region  of  interest 
to  detect  the  tips  of  the  rotor  blades.  The  distance  from  the  rotor  center 
of  rotation  to  the  detected  rotor  tips  provide  the  measurements  necessary 
for  the  final  determination  of  the  rotor  orientation.  If  two  rotor  blade  tip 
locations  can  be  measured,  and  given  the  rotor  center  of  revolution,  and  the 
scaled  length  of  a  rotor  blade,  the  imaged  two-dimensional  elliptical  path 
of  the  blade  tips  can  be  defined.  Given,  this  information  the  tilt  of  the 
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circular  disk  defining  the  rotor  blade  orientation  relative  to  the  helicopter 
fuselage  can  be  determined.  The  two  blade  measurements  required  from 
the  same  image  can  be  either  the  tips  of  two  different  blades  or  the  leading 
and  trailing  edge  of  a  single  blade  during  the  field  integration  time. 


7  Coordinate  Transformations 

In  the  following  sections  the  transforms  are  developed  that  are  necessary  to 
obtain  the  aircraft  orientation  (Euler  angles)  given  the  observer  platform 
and  camera  orientation.  This  discussion  begins  with  the  development  of 
the  basic  transformation  matrices,  then  develops  the  form  of  the  matrix  for 
an  arbitrary  observer  platform  orientation.  The  basic  theory  of  geometric 
transformation  matrices  is  taken  from  Paul  [I4j .  The  treatment  for  a  fixed 
earth  centered  observer  was  taken  from  Brown  [  1  ] .  This  is  then  extended  to 
the  more  general  condition  of  an  arbitrary  (but  known)  observer  orientation 
as  outlined  by  Gonzalez  and  Wintz  [5]. 


7.1  Translation  Transformation  Matrices 

Using  homogenous  coordinates,  the  translation  of  a  vector,  object,  or  coor¬ 
dinate  frame  can  be  performed  using  matrix  multiplication.  The  translation 
to  the  location  described  by  the  vector  w  =  oi  +  6j  +  ck,  is  the  displace¬ 
ment  a,  b.  c  along  the  x,  y,  and  z  axis,  respectively.  Multiplication  by  the 
following  matrix  performs  this  translation: 


Trans{(a,  b,  c))  =  7U-  = 


10  0a 
0  10  6 
0  0  1c 
0  0  0  1 


(5) 


7.2  Rotation  Transformation  Matrices 

Rotations  of  a  point,  object,  or  coordinate  frame  around  a  coordinate  axis 
can  also  be  performed  by  multiplying  the  homogenous  coordinates  by  a 
matrix.  The  rotation  transformation  matrices  can  be  defined  for  rotation 
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about  the  x,  y,  and  7-axis.  A  rotation  through  a  positive  angle  has  the 
convention  that  the  point  is  rotated  counter-clockwise  when  observed  from 
the  point  of  view  of  an  observer  looking  along  the  axis  of  rotation  towards 
the  origin.  The  rotation  matrices  for  rotations  about  the  x,  y,  and  z  axes 
are  shown  below: 


ROT(xJ)  = 


ROT{y,0)  - 


ROT{z,6 )  = 


1  0 

0 

0 

0  cos  8 

- 

sin  0 

0 

0  sin  0 

cos  0 

0 

0  0 

0 
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cos  0 
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sin  6 

0  ' 
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-  sin  8 
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cos  8 
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sin 

0  0 
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cos  0  0 

0 

0 

0  1 

0 

0 

0  0 

1 
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A  sequence  of  rotations  can  be  performed  by  simply  multiplying  the 
appropriate  matrices  together.  However,  the  order  in  when  the  matrices 
are  multiplied  depends  on  on  whether  the  rotation  is  performed  relative  to 
the  base  coordinate  reference  or  the  transformed  coordinates.  A  rotation 
about  the  x-axis  followed  by  a  rotation  about  the  v-axi-  where  this  lat¬ 
ter  rotation  is  about  the  reference  coordinates  is  described  by  the  matrix 
product  ROT(y,6v)ROT(x,Ox).  If  however,  the  rotation  is  about  the  x- 
axis  followed  by  a  rotation  about  (the  now  rotated)  body  y-axis  the  matrix 
product  is  ROT(x.6z)ROT(y,Oif).  The  first  case  describes  the  operations 
performed  in  generating  the  rotated  shape  libraries  for  in  the  image  pro¬ 
cessing  frame  of  reference.  The  later  describes  the  convention  for  the  Euler 
angle  rotations  of  an  airbourne  platform  about  its  body  axes. 


7.3  Reference  Frames 

The  world  coordinate  system  defined  for  aircraft  orientation  is  as  shown 
in  Figure  27.  The  positive  x-axis  points  north,  the  positive  y-axis  points 
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east,  and  the  positive  z-axis  points  down  into  the  plane  of  the  flat  earth. 
The  body  axes  for  the  aircraft  are  embedded  in  the  aircraft  such  that  the 
origin  of  the  body  axes  corresponds  to  the  aircraft  center  of  gravity.  The 
orientation  of  the  aircraft  body  axes  are  defined  as  shown  in  Figure  28. 
The  positive  x-axis  points  out  the  front  of  the  aircraft  nose,  the  positive 
y-axis  point  out  to  starboard,  and  the  positive  z-axis  points  down  through 
bottom  of  the  aircraft. 

The  aircraft  is  in  reference  orientation  when  the  body  axes  and  the  world 
coordinate  axes  are  parallel.  This  corresponds  to  the  aircraft  in  level  flight 
proceeding  north. 

The  coordinate  frame  for  the  image  processing  has  been  defined  differ¬ 
ently.  The  earth  centered  coordinate  axis  for  the  image  processing  is  shown 
in  Figure  29.  In  this  case  the  positive  z-axis  points  north,  the  positive  x- 
axis  points  to  the  west,  and  the  positive  y-axis  points  out  perpendicular  to 
the  flat  earth. 

The  body  centered  coordinates  for  the  image  processing  are  as  showm  in 
Figure  30.  In  this  case,  the  positive  z-axis  points  out  the  nose,  the  positive 
x-axis  points  to  port,  and  the  positive  y-axis  points  out  the  top  of  the 
aircraft. 

When  the  aircraft  is  in  level  flight,  proceeding  north,  the  earth  centered 
coordinate  axes  and  the  body  centered  coordinate  axes  are  parallel. 

In  order  to  describe  points,  objects,  and  coordinate  frames  in  the  world 
coordinate  system,  it  will  be  necessary  to  transform  the  image  processing 
coordinates.  A  rotation  of  the  world  coordinates  by  y-  about  the  x-axis 
followed  by  a  rotation  about  the  original  z-axis  by  y  will  transform  the 
world  coordinate  system  into  the  image  processing  coordinates.  This  trans¬ 
formation  can  be  described  by  the  matrix  Te  where 


Ti:  =  ROT{z.~)ROT{  i,y) 


0  0  10 

-1  0  0  0 

0-100 
0  0  0  1 


(9) 


Premultiplying  a  point,  object,  or  coordinate  frame  by  the  matrix  Te  trans¬ 
forms  the  image  processing  coordinates  to  the  world  coordinates. 
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Figure  29:  Image  processing  coordinate  reference  system 
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7.4  Imaging  Transformations 

Assuming  that  the  camera  is  centered  at  the  origin  of  the  world  coordinate 
reference  system  and  that  the  body  axes  are  parallel.  The  library  of  views 
of  the  aircraft  are  formed  by  rotating  the  object  points  about  the  x,  y,  and 
z  image  processing  coordinate  axes.  Let  Tz ,  Ty.  and  Tz  be  the  rotation 
transformations  about  the  x,  y,  and  z  axes,  respectively.  That  is  Tz  — 
ROT(x,9x ),  Tv  =  ROT{y,$y),  and  Tz  —  ROT (r,  —0 .) .  The  rotation  about 
the  z  axis  is  —  6Z,  since  the  shape  matching  algorithm  measures  the  rotation 
from  the  library  reference  to  the  unknown  for  the  z  axis  rotation. 

The  4  x  P  array  of  the  P  homogeneous  object  coordinates,  O.  is  trans¬ 
formed  as  follows: 

O]  =  TzTyTzO,  (10) 

The  object  coordinates  in  the  image  processing  reference  frame  can  be 
transformed  to  the  world  coordinates  by  multiplying  by  TE  . 

Oi  =  To  O  i 

To  =  TETzTyTt.  (11) 

Finally,  the  object  coordinates  are  projected  on  to  the  image  focal  plane. 
The  final  orientation  of  the  object  reference  frame  is  described  in  the  trans¬ 
formation  matrix  by  the  elements  of  the  matrix.  The  composite  of  all  the 
object  transformations  can  be  described  by  a  single  matrix  having  the  form 

Lz  Mz  Nz  a 
_  Ly  My  A  y  b 

'  Lt  Mz  yz  c 
0  0  0  1 

where  {LZ,LV.  Lz)  is  the  transformed  unit  vector  for  the  object  reference 
x-axis.  Similarly,  (MZ,MV,MZ)  and  (Nz,  Nz)  are  the  transformed  y  and 
z  axes,  respectively.  The  triplet  (a,6.c)  is  the  translation  vector.  The  con¬ 
struction  of  an  homogeneous  transformation  matrix  for  object  coordinate 
frame  orientation  and  position  is  depicted  in  Figure  31. 

The  shape  matching  algorithm  provides  an  estimate  of  the  orientation  of 
the  aircraft  having  undergone  the  object  rotation  transformations.  The  an¬ 
gles  6Z.  8y ,  and  6Z  can  then  be  used  to  determine  the  aircraft  body  centered 
orientation  (i.e.  the  Euler  angle..) 
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Figure  31:  Construction  of  a  homogeneous  transformation 
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7.5  Euler  Angles 

The  orientation  of  an  aircraft  in  the  world  coordinate  system  is  described 
by  its  Euler  angles.  The  Euler  angles  are  the  rotations  of  the  aircraft  about 
its  body  axes.  The  sequence  of  operations  used  to  describe  the  orientation 
is  first  to  rotate  about  the  z-axis  through  the  yaw  angle,  v-  Then  the 
aircraft  is  rotated  about  its  transformed  body  y-axis  through  the  pitch 
angle,  9.  Finally,  the  aircraft  is  rotated  about  its  transformed  body  x- 
axis  through  the  roll  angle,  <f>.  Since  these  rotations  are  performed  on 
the  transformed  coordinates,  the  composite  rotation  transformation  is  the 
matrix  T  —  where  Tv  —  ROT (z ,  v)  ,7*  =  ROT(y,6),  and  T, *.  = 

ROT(x ,  <p). 

The  aircraft  being  tracked  by  the  image  processing  system  has  first  been 
rotated  through  some  unknown  Euler  angles  and  then  the  body  centered 
axes  are  translated  to  the  position  w,.  Therefore,  the  tracked  object  points 
relative  to  the  world  coordinate  system,  Ow,  are  transformed  to 

=  TWtTVtT6,T^,Ow.  (13) 

The  Euler  angles  for  the  tracked  aircraft  can  be  determined  by  equat¬ 
ing  the  elements  of  the  orientation  sub-matrix  for  the  transformation  on 
the  tracked  object,  7’w,7’v.,T'<;1Tf,1  and  the  orientation  sub-matrix  for  the  ob¬ 
served  transformation,  TeTzT^Tz.  This  assumes  the  camera  is  at  the  origin 
of  the  world  reference  system  with  it  coordinate  axes  parallel  with  the  world 
coordinate  axes. 

7.6  Observer  Coordinate  Transformations 

In  most  real  circumstances  the  camera  is  not  located  at  the  world  coordi¬ 
nate  system  origin.  Neither  is  the  camera  coordinate  axes  and  the  world 
coordinate  axes  parallel.  The  camera  is  positioned  onboard  another  air- 
bourne  platform.  As  this  platform  moves  the  camera  orientation  changes. 
It  is  therefore  necessary  to  compensate  for  the  fact  the  the  camera  is  not 
in  the  reference  position. 

First,  the  operations  performed  to  place  the  camera  in  its  actual  position 
and  orientation  will  be  discussed.  Then  the  transformations  necessary  to 
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compensate  for  this  change  observed  orientation  of  the  tracked  object  will 
be  presented. 

The  imaging  system  is  mounted  on  a  gimbal  assembly  at  some  position 
on  the  outside  of  the  aircraft.  Let  G  =  (<7x,  <?,,  <?2)  represent  the  vector 
displacement  of  the  gimbal  assembly  pivot  point  from  the  aircraft  center  of 
gravity.  Let  the  displacement  of  the  imaging  focal  plane  center  from  the 
gimbal  pivot  point  be  described  by  the  vector  f  =  This  trans¬ 

lation  transformation  can  each  be  performed  using  the  following  matrices: 

‘  1  0  0  gz 

Tg=  0  1  0  s* 

0  0  1  gz 

0  0  0  1 

and 

1  0  0  /* 

T  =  0  10/, 

f  0  0  1/, 

0  0  0  1 

The  girnbal  assembly  allows  the  camera  to  be  panned  in  azimuth.  This 
corresponds  to  a  rotation  about  the  z-axis  through  an  angle  6az.  The  camera 
can  be  rotated  in  elevation  corresponding  to  a  rotation  of  angle  6ei  about 
the  y-axis.  Let  Tsa,  and  T represent  the  rotation  transformations  for  the 
rotation  of  the  camera  in  azimuth  and  elevation,  respectively. 

The  aircraft  platform  on  which  the  camera  is  mounted  under  goes  rota¬ 
tion  of  its  body  axes  as  described  by  the  Euler  angle  transformation  ma¬ 
trices  T^cTf,cT^c.  Lastly,  the  position  of  the  platform  relative  to  the  world 
coordinate  system  is  described  as  a  translation  transformation  TWc. 

The  camera  coordinate  frame  therefore  undergoes  the  following  trans¬ 
formation: 

D'c  =  T  0: 

where 

Tc  =  TWc  7V  Tfc  T*r  Tc  Tetl  T6aiTT.  (1C) 

The  transformations  on  the  camera  coordinate  frame  does  not  move  any 
of  the  points  in  the  world  coordinate  system  including  the  points  of  the 
object  being  tracked. 


The  discussion  for  determing  the  tracked  object  Euler  angles  from  the 
observed  rotation  angles  was  based  on  the  assumption  that  the  camera  was 
positioned  at  the  origin  of  the  world  coordinate  system  with  the  axes  paral¬ 
lel  with  the  world  coordinate  axes.  The  camera  can  be  placed  in  standard 
position  by  performing  the  inverse  transformation  on  the  tracked  object 
coordinate  frame.  Therefore,  the  tracked  object  transformation  matrix  is 

Tt  =  Te-1T„tT'iT,tT„.  (17) 

To  determine  the  Euler  angles  given  the  observed  imaging  system  ro¬ 
tation  angles,  the  camera  elevation  and  azimuth  angles,  the  observation 
aircraft  platform  position  &:  orientation,  and  the  displacements  to  the  gim- 
bal  pivot  point  <Lr  gimbal  pivot  to  image  plane,  it  is  necessary  to  equate 
the  elements  of  the  orientation  sub-matrix,  Tt  and  that  of  the  observed 
orientation  sub-matrix,  To- 

Tt  —  1  rp  rp  rp  rp  __ _  rp  rp  rp  rp 

c  -*  vv  ■*  £  -^  z  y  **  i  • 

or 

rv,n,re,  =  rw;r,r£r,rvr„  (is) 

The  elements  of  the  orientation  sub-matrix  of  the  left-hand  side  of  the 
above  equation  are  given  in  Appendix  A.  The  expressions  were  computed 
using  the  symbolic  mathematics  processing  program  REDUCE  [7  . 

Letting  Tt  represent  the  matrix  on  the  right-hand  side  of  the  equation, 
we  can  now  solve  for  the  tracked  objects  Euler  angles. 

So,  we  have 

Ts,T(,T„  =  Tt 

or 

Tf.T,,,,  —  T~'T(.  (19) 

Performing  the  matrix  multiplications,  we  then  have 

cos  Ct  sin  <*>t  sin  9t  cos  c>,  sin  P: 

0  cos  <t>t  —  sin 

-  sin  9t  cos0,  sin<P(  cos  Ct  cos  6t 

Ly  sin  vh  AC  cost,,  -  My  sin  c  ,  .V7  cos  v,  -  A\  s:n  C, 

Lz  sin  tA,  My  cos  v,  -+•  AC  sin  e  ,  Nv  cos  v ,  4-  Nz  sin  v,  (20) 

L*  AC  A'x 
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Now,  equate  corresponding  elements  in  the  two  matrices.  First,  we  have 

Ly  cos  v\  +  Lz  sin  x};t  =  0. 


.  ,i.2  -rLy 

ta  n0,  =  -=- 

Therefore,  we  have  two  possible  solutions  for 

V,'2  =  atan2(TL„,  ±LZ). 
where  v]  and  tf  differ  by  z.  Having, 

sin#(  —  —Lz  and  cos  6t  =  Lz  cos  Vt  —  sin  v’i 

gives 


6t  =  atan2(=fl2,  zz(Z/j  cos  v<1,2  -  sin  y,1'*). 


Lastly, 


(21) 


(22) 


sin  <?*  =  —  ALcosLl  —  ATiSini/>t  and  cos  c?t  =  A 


Therefore, 

Ot  =  atan2(- A”v  cos  b'(1,2  T  Arx  sin  V’,1'2,  :rA/„  cos  Vy'*  ±  A/z  sin  U'/'2)  (23) 

Since  and  d>(  are  dependant  upon  ^>(,  choose  Vi,  and  c><  such  that 

7T  7T 

-z  <  xpt  <  z  and - <  ot  <  —  and  —  7r  <  d>(  <  z. 

2  2 

When  the  aircraft  is  in  vertical  dive  or  climb,  then  Lv  —  Lz  =  0  and  psit 
is  undefined.  In  this  case,  L't  =  0  by  convention. 


8  Conclusions 

In  this  study  the  results  are  presented  for  the  comparing  the  use  of  two 
global  shape  recognition  algorithms  applied  to  the  problem  of  helicopter 
type  classification  and  orientation  determination.  Under  the  conditions  of 
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varying  library  density  and  simulated  imaging  noise,  it  is  shown  that  the 
Fourier  descriptor  method  is  superior  in  performance  measured  as  classifi¬ 
cation  accuracy  and  mean  (or  median)  angle  error. 

In  addition,  a  scheme  is  proposed  describing  outlining  the  library,  seg¬ 
mentation,  and  processes  for  *he  determination  of  rotor  tip  path.  The  rela¬ 
tive  orientation  of  the  blade  path  disk  can  be  determined  from  the  elliptical 
shape  of  the  projected  rotor  tip  path. 

Finally,  the  transformations  for  computing  the  tracked  object  Euler  (ori¬ 
entation)  angles  are  derived. 


61 


9  References 


1.  Robin  E.  Brown,  “Transformation  of  Observed  Aircraft  Orientation  to 
Euler  Angles,”  l'.S.  Army  Summer  Assoc.  Program  for  High  School 
Science  and  Mathematics  Faculty,  Aug.  22,  198G. 

2.  Charles  R.  Giardina  and  Frank  P.  Kuhl,  “Accuracy  of  Curve  Approx¬ 
imation  by  Harmonically  Related  Vectors  with  Elliptical  Loci.”  Com¬ 
puter  Graphics  and  Image  Processing  6,  pp.  277-285,  1977. 

3.  Charles  R.  Giardina,  Frank  P.  Kuhl,  T.A.  Grogan,  and  O.  Robert 
Mitchell,  “A  Spatial  Domain  Walsh  Feature  Set,”  Trans,  of  the  28'”' 
Conference  of  Army  Mathematicians ,  ARO  Report  83-1,  pp.  397-408, 
February,  1983. 

4.  Marcus  Elgin  Glenn,  “Fourier  Descriptors  of  Two  and  Three  Dimen¬ 
sional  Objects,”  M.S.E.E.  Thesis,  Purdue  University,  May,  1982. 

5.  Rafael  C.  Gonzalez  and  Paul  Wintz,  Digital  Image  Processing.  Second 
Edition,  Addison-Wesley,  1987. 

G.  Timothy  A.  Grogan,  “Shape  Recognition  and  Description:  A  Compar¬ 
ative  Study,”  Ph.D.E.E.  Thesis.  Purdue  University,  August,  1983. 

7.  Anthony  C.  Hearn.  “REDUCE  User's  Manual,  Version  3.0",  Prim! 
Publication  CP78(-l/83),  19S3. 

8.  P.V.C.  Hough,  “Methods  and  Means  for  Recognizing  Complex  Pat¬ 
terns,”  U.S.  Patent  3,0G9,G54,  19G2. 

9.  Jane's  All  the  Worlds  Aircrafi.pp.  335,  337,  4G2,  4G3,  510,  51 1,  (1985). 

10.  Frank  P.  Kuhl  and  Charles  R.  Giardina,  “Elliptic  Fourier  Features  of  a 
Closed  Contour,”  Computer  Graphics  and  Image  Processit.g,  \oi.  18, 
pp.  230-258,  1982. 

11.  Mac3D,  Version  2.0,  Challenger  Software  Corporation.  1985. 

12.  O.  Robert  Mitchell.  Frank  P.  Kuhl.  Timothy  A.  Grogan,  and  Didur 
Charpcnticr,  “A  Shape  Extraction  and  Recognition  System,”  South- 
con/SL’ ,  March  23-25,  Orlando,  Florida.  pp.(  l  1)  1-1.  19v2. 


13.  O.  Robert  Mitchell.  Anthony  P.  Reeves,  and  Timothy  A.  Grogan,  “Al¬ 
gorithms  and  Architectures  for  Global  Shape  Analysis  in  Time- Varying 
Imagery,”  SP IE  Proceedings  SCO,  Robotics  and  Industrial  Applications , 
San  Diego,  Aug.  24-27,  1982. 

14.  Richard  P.  Paul,  Robot  Manipulators:  Mathematics,  Programming, 
and  Control,  MIT  Press,  1981. 

15.  E.  Persoon  and  K.S.  Fu,  “Shape  discrimination  using  Fourier  descrip¬ 
tors,”  IEEE  Trans.  Syst.,  Man,  Cybcrn.,  Yol.  SMC-7,  pp.  170-179, 
Marc!:.  1977. 

16.  Timothy  P.  Wallace  and  P.A.  Wintz,  “An  efficient  three-dimensional 
aircraft  recognition  algorithm  using  normalized  Fourier  descr  ptors,” 
Computer  Graphics  and  Image  Processing ,  Vo).  13,  pp.  99-126,  1980. 

17.  Timothy  P.  Wallace  and  Owen  R.  Mitchell,  “Analysis  of  three-dimensional 
movement  using  Fourier  descriptors,”  IEEE  Trans.  Pattern  Anal,  and 
Mach.  Intcll.,  Yol.  PAM1-2,  no.  6,  pp.  583-588.  November,  1980. 


64 


APPENDIX  A 

Raw  Data  for  Classification  and  Orientation  Accuracy 
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The  following  pages  contain  the  raw  data  for  the  experiments  determining 
the  classification  and  orientation  accuracy  for  the  helicopter  fuselages.  The 
information  is  organized  into  three  groups  results:  those  for  the  Fourier 
descriptors,  the  Walsh  points,  and  Fourier  descriptors  using  multiple  nor¬ 
malizations.  Then  for  each  group,  the  results  are  tabulated  under  the 
variations  in  the  library  density  and  the  “unknown”  image  signal-to-noise 
rat  10. 

The  libraries  are  numbered  1,  2,  and  3,  corresponding  to  the  Bell  Model 
203  ril-1.  McDonnell!  Douglas  500MD,  and  Sikorsky  LTIGOA  Black  Hawk, 
respectively. 
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l.ibraiy  Size  <)xJ  SNR  3  Features  Fourier  rnoffs 


Confusion  Matrix 
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42*  <-<»rro<‘t  i'|v.r.i!ic,\l  ion  ISM/  ]  MM  vhirh  is 


Confusion  Matrix 
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Library  Size  9x/  SNR  3  Features  Walsh  Points 


t  ho 
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t’  l  .)  4  iV.»A  1 1)  ftUl  q.H*. 
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Library  Siz«  lx  >  sw  »  K.-ii  ur<"-.  Fourier  (Mult.  N-.iml/.n  J  Library  Si?.<-  ■>*  1  SNF  1  (Vat  urr»;  F<mri»r  (Mult.  N'umlrn 


APPENDIX  B 


Object  Orientation  Matrix 
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Below  are  the  expressions  for  the  elements  of  the  object  orientation  matrix 
as  computed  using  the  REDUCE  [7]  program.  The  orientation  matrix  is 
the  upper  3x3  submatrix  of  the  matrix  formed  by  the  transformation 


J*  rT~'  'T*  ep  rj~y  p  nn  <p  rr-  ryi  rp  rp 

w(JwcUc-'«d^Jch,ih,jJrJ£:hhh 


The  correspondence  between  the  parameters  of  the  above  transformation 
and  the  symbols  used  in  the  listing  are  shown  below: 

0Z  =•  thetax 

0V  =  thetay 

6Z  —  thetaz 

ty  =  psic 

0,  —  thetac 

<p:  =  phic 

OJ  =  thetael 
6a:  =  thetaaz 


Lx  =  cos (phic) *cos (psic) ‘cos (thetaaz) *  Bin (thetac) *  sin (thetael) * 
Bin (thetay)  +  cos (phic) ‘cos (psic) ‘cos (thetael) ‘cos (thetay) 
‘sin(tbetae) ‘sin(thetaz)  +  cob (phic) *cos (psic) *cob (thetay) 
*cos(thetaz)*8in(thetaaz)*Bin(thetac)*8in(thetael)  -  cos( 
phic) ‘cob (thetaaz) ‘cos (thetay) ‘cos (thetaz) ‘sin (psic)  +  cos 
(phic) *8in(paic) ‘sin (thetaaz) *sin (thetay)  -  cos (psic) *cos ( 
thetaaz) »cos (thetac) ‘cos (thetael )» sin (thetay )  -  cos(p8ic)‘ 
co 8 (thetaaz) * cos (thetay) ‘cos (thetaz) *  sin (phic) ‘Bin (thetac) 
-  cos (psic) * cos (thetac) ‘cos (thetael) ‘cos (thetay) * cos ( 
thetaz) ‘Bin (thetaaz)  +  cob (psic) ‘cos (thetac) ‘cob (thetay) • 
sin (thetael )* b in  (thetaz)  +  cos (psic )» s in (phi c )* s in (thetaaz 
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) ‘sin(thetac) ‘sin(thetay)  -  cob (thetaaz) *ain (phic) »ein( 
psic) ‘sin(thetael) ‘ain(thetay)  -  cob (thetael) ‘cob (thetay ) * 
ain(phic) ‘sin(psic) ‘Bin(thetaz)  -  cob (thetay) *cob (thetaz) * 
sin (phic) * Bin (paic) *sin(thetaaz) *  Bin (thetael) 

-  cos (phic) *cob (psic) *cob( theta az)» cob (t he t ay )*ein(thetac 
) ‘sin(thetael) ‘sin(thetax)  -  cos (phic) ‘cos (psic) ‘cos  ( 
thetael) ‘cob (thetax) ‘cos (thetaz) ‘Bin (thetac)  +  cos(phic)» 
cob (psic) *c os (thetael) *ain (thetac ) * sin (thetax) *  Bin (thetay) 
‘ain(thetaz)  ♦  cob (phic) *cos (psic) ‘cos (thetax) ‘sin (thetaaz 
)*sin(thetac)*Bin(thctael)*sin(thetaz)  +  cob (phic) ‘cob ( 
psic) * cob (thetaz) * sin (thetaaz) ‘Bin (thetac) ‘sin (thetael)* 
Bin(thetax) ‘Bin(thetay)  -  cob (phic) ‘cob (thetaaz) ‘cos ( 
thetax) ‘Bin(psic) ‘ain(thetaz)  -  cob (phic) *cob (thetaaz) ‘cos 
(thetaz) »sin (psic) ‘Bin (thetax) ‘Bin(thetay)  -  cos(phic) ‘cos 
(thetay)*8in(p8ic)«Bin(thetaaz)*Bin(thetax)  +  cos(pBic)* 
cos (thetaaz) *cob (thetac) * cob (thetael) *coa (thetay) ‘sin ( 
thetax)  -  cos (psic) ‘cob (thetaaz) ‘cos (thetax) ‘Bin (phic) ‘Bin 
(thetac)*sin(thetaz)  -  cob (psic) ‘cob (thetaaz) *cob (thetaz) » 
Bin(phic)‘8in(thetac)»Bin(thetax)»sin(thetay)  -  cos(pBic)» 
cos (thetac) ‘cos (thetael) *cos (thetax) ‘sin (thetaaz) ‘Bin ( 
thetaz)  -  cob (psic) ‘cos (thetac) ‘cos (thetael) ‘cob (thetaz) * 
8in(thetaaz)*Bin(thetax)*sin(thetay)  -  cob (psic) ‘cos ( 
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thetac) *cos (thetax) *cos (thetaz) *  sin (thetael)  +  cos(psic)* 
cos(thetac) *  sin (the tael) *  sin (thetax) *sin(thetay)*sin( 
thetaz)  -  cos (psic) *cos (thetay )* sin (phic ) *sin (thetaaz) *sin 
(thetac) *sin(thetax)  +  cos (thetaaz) *cos (thetay) *sin (phic ) * 
sin  (psic)  *  sin  (thetael)  *sin  ''thetax)  +  cos (thetael ) *cos  \ 
thetax) *cos (thetaz) *8 m (phic) *sin (psic)  -  cos (thetael )* sin 
(phic) *sin (psic) *sm (thetax) *  sin (thetay) *  sin (thetaz)  -  cos 
(thetax) »sin (phic)* sin (psic) • sin (thetaaz)* sin (thetael )*sin 
(thetaz)  -  cos (thetaz) *sin (phic ) »sin (psic) *  sin (thetaaz) » 
sin (thetael) * sin (the tax) *  sin (thetay) 

-  cos (phic )*cos (psic) *cc  8 (thetaaz) *cos  (thetax) *cos (thetay 
) *sin(thetac) *sin (thetael )  +  cob (phic) *cos (psic) *cos ( 
thetael )*cos (thetax) « sin (thetac) *sin (thetay)* sin (thetaz) 

+  cos (phic) *cos (psic) *cos (thetael) *cos (thetaz) *sin (thetac 
) *sin (thetax)  +  co- (phic) *cos (pBic ) *cos (thetax) *cos (thetaz 
) *6in (thetaaz) *sin (thetac) *sin (thetael) *sin (thetay)  -  cost 
ph ic)*cos (psic) *s in (thetaaz)* sin (thetac)* sin (thetael)*sin( 
thetax)  *sm(thetaz)  -  cos  (phic)  »cos  (thetaaz)  *cos  (thetax)  * 
cos (thetaz) * s in (psic) *  sin (thetay )  +  cos (phic) *cos (thetaaz) 
*sin(p3ic) *sin(thetax) *Bin(thetaz)  -  cos (phic ) *cos (thetax) 
•  cos (thetay )* sin (psic )* sin (thetaaz)  +  cos (psic) *cos  ( 
thetaaz )* cos (thetac )*cos (thetael)* cos (thetax)* cos (thetay) 
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-  cos  (ps  ic  )  *  coa  ( theta  az  )  *cos  (the  tax)  *cos  (thetaz)  *  sin  (pkic 
) *sin(thetac) *sin(thetay)  +  cos (psic) *cos (thetaaz) *Bin ( 
phic) *sin(thetac) *sin(thetax) *sin(thetaz)  -  cos (psic) *cos ( 
thetac) *cos (thetael ) *cos(thetax)*cos(thetaz)*sin (thetaaz) * 
sin(thetay)  +  cos (psic) *cos (thetac) *cos (thetael) *sin ( 
thetaaz) *sin(thetax) *sin(thetaz)  +  cos (psic) *cob (thetac ) » 
cos (thetax) *Bin (thetael) *sin (thetay) *sin(thetaz)  +  cob( 
pflic)*co8 (thetac) *coa(thetB2) *ein (thetael) *ain (thetax)  - 
cob (psic )» cos (thetax) *cos (thetay)* sin (phic) *  sin (thetaaz) * 
sin (thetac)  +  cos (thetaaz) *cos (thetax) *cos (thetay) *sin( 
phic) *sin (psic) *sin (thetael)  -  cos (thetael) *cos (thetax) * 
sin(phic)*Bin(psic)*Bin(thetay)*Bin(thetBz)  -  coB(thetael) 
•coa (thetaz) *ein (phic ) *ein (peic) *ein (thetax)  -  coe(thetax) 
•cos (thetaz) « sin (phic)*e in (psic)* sin (thetaaz) • sin (thetael) 
•sm(thetay)  +  sin  (phic)  •  s  in  (ps  i  c  )*  sin  (thetaaz)  *  s  in  ( 
thetael)* sin  (thetax)* sin  (thetaz) 

-  cos (phic) *cos (psic) *cos (thetaaz) *cos (thetay) *cos  (thetaz 
)  +  cos (phic) *cob (psic) *sin (thetaaz) *  sin (thetay )  -  cost 
phic)*cos(thetaBz)*Bin(psic)*Bin(thetac)*sin(thetael)*Bin( 
thetay)  -  cos (phic) *cos (thetael ) ‘cos (thetay )* sin (psic )* sin 
(thetac) *Bin(thetaz)  -  cos (phi c) *cos (thetay ) *cos (thetaz) * 
sin(psic) » sin (thetaaz) *sin(thetac) *sin(thetael)  -  cos(pBic 
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)  *  cos (thetaaz)  *  sin (phic) * s in (thet ael ) • sic (thetay )  -  cost 
psic)  *cos (thet ael ) *cos (thet ay ) *  Bin (phic) * ein (thetaz)  -  cob 
(ps ic )*cos (thetay )*cos(thetaz)* sin (phic)*8in (thetaaz)* sin ( 
thetael)  +  cob ( the t aaz) *  cob (thetac ) • cob (thet ae 1 ) * s in  (psi c ) 
•sin (thetay)  +  cob (thetaaz) *cob (thetay ) *cos (thetaz) *ain ( 
phic)  *8in(psic)  *sin(thetac)  +  cos  ( the  t  ac )  »c  os  ( thet  ae  l  'i  *  c  os 
(thetay) *cos (thetaz) *ain (psic) •• in (thet aaz)  -  cos(thetac)* 
cob (thet ay )• e in (peic) *ain (thetael) *  sin (thetaz)  -  ain(phic) 
*ein(psic) *s in (thetaaz) *  sin (thetac) *Bin(thetay) 

-  cos (phic)*cos (psic )*cos (thet aaz) *cos(thetax)*Bin (thetaz 
)  -  cob (ph ic )*coa (psic )*cos (thetaaz)* cob (the taz)*B in ( 
thetax) *sin(thetay)  -  cob (phic) *cos  (psic) »cob  (thetay) *sm  ( 
thetaaz) *sin  (the tax)  +  cos (phic )*cos (thetaaz) «cos (thet ay)* 
■  in (psic ) *sin (thetac) *  sin (thetael) *Bin (thetax)  +  cos(phic) 

•  cos (thetael )*cos (thetax) *cos (thetaz) *sin(pBic)» sin  (thetac 
)  -  cos  (phic) *cos (thet ael) *sin (psic) *  Bin (thetac) *  sin ( 
thetax) *ein(thetay) *sin(thetaz)  -  cos (ph  i  c ) *  cob  ( thet ax ) * 
Bin(psic)*sin(thetaa2)*Bin(thetac)*Bin(thetael)*8in(thetaz 
)  -  ccs(phic)*co8(thetaz)*Bin(pEic)*sin(thetaaz)*8ii;( 
thetac) *sm (thetael) • sin (the tax)*  Bin (thetay)  +  cob (psic)* 
cos(thetaaz)*co8(tbetay)*Bin(phic)*8in(thetael)*Bin(thetax 
)  ♦  coB(psic)*coB(thetael)*cos(thetax)*cos(thetaz)*8in( 


phic)  -  cos (psic) ‘cos (thetael) ‘sin (phic) ‘Bin (thetax) ‘sin ( 


thetay) ‘sin(thetaz)  -  cob (psic) ‘cos (thetax) ‘sin (phic) *ein ( 
thetaaz) »sin (thetael ) ‘ain(thetaz)  -  cos (psic ) ‘cos (thetaz) * 
sin(phic)»sin(tbetaaz)*sin(thetael)»sin (thetax) ‘sin (thetay 
)  -  cos (thetaaz) ‘cos (thetac) »cos (thetael) »cos (thetay) »sin ( 
psic) • s in (thetax)  +  cos (the t aaz) »cos (thetax) » sin (phic )» Bin 
(peic) »sin (thetac) ‘ain(thetaz)  +  cos (thetaaz) ‘cos (thetaz) » 
Bin (phic) ‘sin (psic) ‘Bin (thetac) ‘Bin (thetax) *  Bin (thetay)  + 
cos (thetac) ‘cos (the tael) ‘cos (thetax) ‘sin (psic) ‘Bin (thetaaz 
) «ein (thetaz)  +  cos (thetac) ‘cos (thetael) »cos (thetaz) ‘sin ( 
psic) ‘sin (thetaaz) »sin (thetax) ‘sin (thetay)  +  cos(thetac)» 
cos(thetax)‘cos(thetaz)*sin(p8ic)»sin(thetael)  -  cos( 
thetac) ‘sin (psic) *sin (thetael) ‘sin (thetax) ‘sin (thetay) »sin 
(thetaz)  +  cob (thetay) ‘sin (phic) ‘sin (psic) ‘sin (thetaaz) * 
sin( thetac) ‘sin(thetax) 

-  cos (phic) ‘cos (psic) ‘cos (thetaaz) ‘cos (thetax) ‘cos (thetaz 
) ‘sin (thetay)  +  cos (phic) ‘cos (psic) ‘cos (thetaaz) ‘sin ( 
thetax) ‘sin(thetaz)  -  cos (pfc ic) ‘cos (psic) ‘cos (thetax) ‘cos ( 
thetay) ‘sin (thetaaz)  +  cos (phic) ‘cos (the taaz)*cos (thetax)* 
cos (thetay) »sin (psic) ‘sin (thetac) ‘sin(thetael)  -  cos(phic) 
*coB(thetael)»cos(thetax)*sin(psic)‘sin(thetac)*8in(thetay 
) ‘sin (thetaz)  -  cos (phic ) ‘cos (thetae 1 ) ‘cos (thetaz) *  sin ( 
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psic) »sin (thetac ) ‘sin  (thet ax)  -  cob (phic) ‘cos (thetax) ‘cob ( 
thetaz) *sin(psic)»Bin(thetaaz)*Bin (thetac) *ein (the tael) * 
sin(thetay)  +  cos (phic )* s in (ps ic) *  Bin (thetaaz) * s in (thetac ) 
*sin (thetael) ‘sin (thetax) ‘Bin (thet az)  +  cos (psic) »cob ( 
thetaaz) *cos (thetax) »cos(thetay) *sin (phic)* sin (the tael)  - 
cos(psic)»cos(thetael)»co6(thetax)»Bin(phic)»Bin(theiay)» 
sin(thetaz)  -  cos (psic) *coa (thetael) ‘cos (thetaz) »ain (phic) 
•sin(thetax)  -  cob (psic) »cob (thetax) »cob (thetaz) »ein (phic) 
♦sin (thetaaz) ‘sin (thet ael) ‘Bin (thetay)  +  cob (psic) »sin ( 
phic)»sin(thetaaz)‘ein(thetael)*s in (thetax) »sin (thetaz)  - 
cos (thetaaz) »cos (t he tac)»cos (thetael) ‘cob (thetax) ‘cob ( 
thetay)‘8in(psic)  +  cob (thetaaz) ‘cob (thetax) ‘cos (thetaz) * 
Bin (phic ) ‘sin (psic) ‘sic (thet ac) *  Bin (thet ay )  -  cos(thetaaz) 
•sin (phic) ‘sin (psic) ‘Bin (thetac) *8in (thetax) ‘Bin (thetaz) 

‘  cos (thetac) *co8(thetael)*coB (thetax)* cob (thetaz) ‘sin ( 
psic) ‘Bin(thetaaz) ‘Bin(thetay)  *  cos (thetac) ‘cob (thetael) * 
sin (ps ic) ‘sin (thetaaz) *ein (thetax) *  sin (thetaz)  -  cos( 
thetac) ‘cos (the tax) ‘sin (psic) »sin (thetael) ‘Bin (thetay )*sin 
(thetaz)  -  cob (the tac) *cob (thet az) » sin (psic ) ‘Bin (thet ael) * 
Bin(thetax)  +  cos ( thet ax ) * c ob ( thet ay ) • s in (phi c ) » s in (ps  i  c ) * 
Bin(thetaaz) ‘Bin(thetac) 
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Lz 


cos(phic)*cos (thetaaz) *cos(thetac) *  Bin (the tael) ‘sin (thetay 


Mz 


)  +  cos (phic) *cos (thetac) *  cob (thetael) *cos (thetay) *sin( 
thetaz)  +  cos (phic) ‘cos (thetac) ‘cos (thetay ) ‘cos (thetaz) » 
sin(thetaaz) *sin (thetael)  -  cos (thetaaz) *cos (thetac) ‘cos  ( 
thetay ) ‘cos (thetaz) ‘sin (phic)  +  cos (thetaaz) ‘cos (thetael) * 
sin(thetac) »sLn(thetay)  +  cos (thetac) *ein (phic) *sin ( 
thetaaz) *sin(thetay)  +  cos (thetael) *cos (thetay) *cos (thetaz 
) *sin (thetaaz) *sin (thetac)  -  cos (thetay ) *sin(thetac) *sin( 
t he tael)*s in (thetaz) 

-  cos(phic)*cos (thetaaz) ‘cos (thetac) * cos (thetay) *Bin ( 
thetael) »sin(thetax)  -  cos (phic) *cob (thetac) ‘cos (thetael ) * 
cos(thetax)»coa(thetaz)  +  cos (phic) ‘cob (thetac) ‘cos ( 
thetael)*sin(thetax)*8in(thetay)»8in(thetaz)  +  cos(phic)» 
cos (thetac) ‘cos (the tax) * sin (thetaaz) ‘sin (thetael) ‘sin ( 
thetaz)  +  cos (phic) ‘cos (thetac) » ros (thetaz) ‘sin (thetaaz) * 
sin (thetael) ‘sin (thetax) ‘sin  (thetay )  -  cos(the  aaz)*cos( 
thetac) ‘cos (thetax) ‘sin (phic) ‘sin (thetaz)  -  cos (thetaaz) » 
cob (thetac) ‘cos (thetaz) ‘sin (phic) ‘sin (thetax) ‘sin (thetay) 

-  cob (t he ta az ) ‘cob (the tael) »cos (thetay) ‘sin (thetac) *s in ( 
thetax)  -  cos (thetac) ‘cos (thetay ) ‘sin (phic )* sin (thetaaz) • 
sin(thetax)  +  cos (thetael) ‘cos (thetax) ‘sin (thetaaz) ‘sin ( 
thetac) ‘sin(thetaz)  +  cos (thetael) ‘cos (thetaz) »sin (thetaaz 
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) *sm (thetac) *sin (thetax) *sin (thetay )  +  cos (thetax) *cos ( 


thetaz) *sin(thetac) *sin(thetael)  -  s in (thetac )* sin (thetael 
)* sin (thetax) *sin (thetay ) »  sin (thetaz) 

!,'z  =  -  cos (phic) *cos (thetaaz) *cos (thetac) *cos (thetax) *cos  ( 

thetay) *sin (thetael)  +  cos (phic )* cos (thet ac) *  cos (thetael) * 
cos (thetax) *sin (thetay ) *sin (thetaz)  +  cos (phic) *cos (thetac 
) *cos (thetael ) *cos (thetaz) *sin (thetax)  +  cos (phic) *cos ( 
thetac)* cos (thetax)* cos (thetaz)* sin (thetaaz)* sin  (thetael)* 
sin(thetay)  -  oos (phic) *cos (thetac) »sin (thetaaz) *sin ( 
thetael) *sin (thetax) »sin(thetaz)  -  cos (thetaaz) *cos (thetac 
) *cos (thetax) *cos (thetaz) *Bin (phic) *sin (thetay)  +  cob( 
thetaaz)* cos (thetac) *s in (phic) *ein ( thetax ) *ain  (thetaz)  - 
cos (thetaaz) ‘cos (thetael) *cos( thetax) *cos (theta y)*sin( 
thetac)  -  cos (thetac) ‘cos (thetax) *cob (thetay) *  sin (phic) * 
sin(thetaaz)  *  cos ( thetae 1 ) *cob ( the t ax) * c os ( thet ez ) * b i n ( 
thetaaz)*sin(thetac)*sin( thetay)  -  coB(tbetael)*sm( 
thetaaz) *  sin (thetac) *  sin (thetax) *  sin (thetaz)  -  cos(thetax) 

* sin(thetac) *  sin (thet ael) *sin (thet ay) *  Bin (thetaz)  -  cost 
thetazl*sin(thetac)*sin(thetael)*sin(thetax) 


or, 
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